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We calculate exactly the two point mass-mass correlations in arbitrary spatial dimensions in the 
5— ( , aggregation model of Takayasu. In this model, masses diffuse on a lattice, coalesce upon contact 

and adsorb unit mass from outside at a constant rate. Our exact calculation of the variance of 
' mass at a given site proves explicitly, without making any assumption of scaling, that the upper 

critical dimension of the model is 2. We also extend our method to calculate the spatio-temporal 
' correlations in a generalized class of models with aggregation, fragmentation and injection which 

include, in particular, the g-model of force fluctuations in bead packs. We present explicit expressions 
for the spatio-temporal force-force correlation function in the g-model. These can be used to test 
' the applicability of the g-model in experiments. 
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I. INTRODUCTION 



. It is well established that interacting many body non-equilibrium systems, evolving via the dynamics of their 
^ microscopic degrees of freedom, can reach a large variety of steady states in the limit of large time. Among these 
I I steady states, those characterized by power law distributions of different physical quantities, over a wide range of 
parameter space, have received a lot of attention in the last decade. The phenomenon of the emergence of power 
law distributed steady states without fine tuning of external parameters has been dubbed self organized criticality 
O ; (SOC) §]. The concept of SOC has been very useful in understanding scale invariance in a large number of physical 
systems including sand piles, driven interfaces, river networks and earthquakes. One of the early models of SOC is 
the mass aggregation model proposed by Takayasu . This simple lattice model describes a system in which masses 
^ diffuse, coalesce upon contact and adsorb unit mass from outside at a constant rate. In the limit of large time, the 
CNj mass distribution evolves into a time- independent form with a power law tail for large mass, in all dimensions. This 
\l , time-independent mass distribution was computed exactly in one dimension and within mean field theory |]. The 
^ ' appearance of a power law distribution without fine tuning of parameters and the fact that the power law exponent 
, can be computed exactly makes the Takayasu model (TM) one of the simplest analytically tractable models of SOC. 
' TM also has close connections |3| to other models of statistical mechanics such as the Scheidegger river network 
, model the directed abelian sandpile model |^ and the voter model TM is exactly equivalent to Scheidegger 's 
river model which had been proposed to explain theoretically the observed drainage patterns of river catchment area, 
in particular Hack's law which relates the catchment area to the length of the river ||^. Though more sophisticated 
models have been proposed to explain Hack's law, Scheidegger's model is a simple stochastic model which explains 
the law with fair accuracy. In the case of the directed abelian sandpile model the probability of having an avalanche 
of size m turns out to be identical to the mass distribution P{m) in the TM. The voter model space time trajectory is 
O ■ very similar to that of the TM except that the processes proceed in the opposite direction in time. Thus, a complete 
^ ' solution of the TM would also help in understanding these related models better. 

^ . There has been a recent revival of interest in the TM as simple modifications of the model led to the understanding 
' of a variety of systems including force fluctuations in granular media such as bead packs [||, various aspects of river 
i^jj networks particle systems in one dimension with long range hops fictl and generalized mass transport models 
d ' PI . In addition, two other natural generalizations of the TM were found to display nontrivial nonequilibrium phase 
transitions in the steady state [p^p^ . 

Even though the single site mass distribution function in the TM has been computed exactly in 1-dimension and 
in the mean field limit, the spatial and temporal correlations between masses have remained an open question. In 
fact, Takayasu and Takayasu, in their recent review article JT^ , have commented on the difficulty of computing the 
spatial mass-mass correlation function both analytically and numerically. In this paper, we calculate exactly both the 
spatial and temporal mass-mass correlation functions in the TM in all dimensions. The calculation of the variance of 
mass in all dimensions also settles rigorously that the upper critical dimension of the TM, beyond which the mean 
field exponents are correct, is 2. 
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Our technique can also be used to calculate the spatio-temporal two-point correlations in a class of models which 
are generalizations of the TM. In particular we calculate exactly the two-point force-force correlation function in the 
g- model of force fluctuations in bead packs If experiments can be devised to measure these spatio-temporal 

correlations in bead packs, then our exact results would be useful for comparison. 

The paper is organized as follows. In section II, we define the TM and briefly review the earlier results on this model. 
We also review the known results for the g-model of force fluctuations in bead packs and then summarize our main 
results. In section III, we study the spatio-temporal correlations in the TM in one dimension. In subsection III-A, we 
compute exactly the equal-time mass- mass correlations between two points in space for all times. We also perform 
a Monte Carlo simulation to compute this correlation function numerically and find excellent agreement between 
numerical and exact results. In subsection III-B, we solve exactly the temporal autocorrelation and derive its large 
time behavior. In section IV, we study the spatio-temporal mass-mass correlations in arbitrary dimensions. In section 
V, we study exactly the spatio-temporal correlations in the generalized g-model. We present the explicit results for 
the correlation functions in the g-model when the distribution of transmitted fractions of weights are uniform. We 
also provide numerical evidence for the equal-time correlations. Finally, we conclude in section VI with a summary 
and outlook. 



For simplicity we define the TM on a one dimensional lattice with periodic boundary conditions. Generalization 
to higher dimensions is straightforward. Each site of the lattice has a nonnegative mass variable. Given a certain 
configuration of masses at time t, the system evolves via the following dynamics. Evolution at each discrete time 
step consists of two moves : (1) with probability 1/2 each mass hops to its right and with probability 1/2 it stays 
at the original site and (2) a unit mass is added to each site. The first move corresponds to the diffusion of masses 
while the second move corresponds to injection of unit masses from outside. If the diffusion move (1) results in two 
masses coming to the same site, then the total mass at the site simply adds up. Thus, the evolution of the masses is 
described by the stochastic equation. 



where the 's are independent and identically distributed random variables taking values or 1 with equal probability 
1/2, and rrii is the mass at site i. While the first move (1) tends to create big masses via diffusion and aggregation, 
the second move (2) replenishes the lower end of the mass spectrum. The competition between the two, leads in the 
long time limit, to a time-independent single-site mass distribution with a power law tail for large mass. This happens 
irrespective of the initial condition. For convenience, one can start from an initial configuration that has zero mass 
at each site. 

Note that the dynamics of the TM defined above is parallel, i.e., all sites are updated simultaneously in every time 
step. Alternately one could define the model in continuous time where the mass at every site hops to the right with 
rate p and injection of unit masses at every site occurs with rate 1. It turns out that the large distance and long time 
behaviors of the TM are insensitive to the particular type of dynamics used. This is in contrast to other recently 
studied generalized mass models [ ^l)|JTl| ] where the steady state mass distribution depends non universally on the type 
of dynamics. It turns out that while the parallel version of the TM is convenient for establishing the mapping to 
other models, the continuous time version is sometimes easier for the purpose of calculations. We will therefore use 
either of the two versions whichever is more convenient in a particular situation. For example, in one dimension we 
will calculate the mass-mass correlation function with parallel dynamics while for arbitrary dimensions d, we will use 
the continuous time version as significant simplifications occur in that case. 

We now briefly review the earlier results on the TM. Takayasu and coworkers [p|jl^] originally showed that the 
probability P{m, t) that a site has mass m at time t approaches a time-independent form in the long time limit 
t oo. This time-independent distribution was shown to have a power law tail, ~ m""^ for large m and the exponent 
T was computed exactly in d = 1 (rid = 4/3) and within mean field theory (rm/ — 3/2). It was also shown that 



II. MODEL AND RESULTS 



A. Takayasu Model 



■mi{t){l - Ti) + mi-i{t)ri_i + 1 



(1) 




(2) 
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where the exponent S is related to r via the simple scaling relation, 5 = 1/(2 — r) ||T^. Recently Swift et. al. 
have argued that in d spatial dimensions, r = 2{d+l)/{d + 2) for d < 2 and r = 3/2 for d > 2. In d = 2, they argued 
that in the limit t ^ oo, P{m) ^ r7i~^/^(logm)^/^ for large m indicating that 2 is the upper critical dimension of the 
TM. It is also known ||] that the power law distribution of mass is stable with respect to fluctuations in the initial 
conditions and is insensitive to whether the particles hop symmetrically in space or not. 

While the single site mass distribution in the TM is known analytically as mentioned above, the spatio-temporal 
correlations between masses are far from being understood. Recently Takayasu and Takayasu have pointed out that 
while there exist spatial correlations in TM, they are difficult to compute not only analytically but even numerically 
. In addition there is no rigorous derivation of the upper critical dimension dc of the model. While equivalence 
with the voter model would predict dc = 2 early numerical simulations suggested dc — 4. The argument of 
Swift et. al. ||T^ that dc = 2 is also not rigorous as it relied on a scaling ansatz for the age distribution of particles 
which looked plausible but was not proved. 

B. The q-Model of Force Fluctuations 

The g-model was proposed [^0as a simple scalar model to understand the distribution of forces observed in 
real three dimensional bead packs |lq] . This model assumed that the force chains observed in experiments were due 
to inhomogeneity in packing leading to unequal distribution of weights supported by a bead. Ignoring the spatial 
correlations between inhomogeneities, the model considered a regular lattice of sites, each containing a bead of mass 
unity. The total weight on a given bead at a layer is transmitted randomly to 2 nearby beads in the layer underneath 
(In the original version of the model, N adjacent beads were considered). Let m{i,t) be the weight supported by a 
bead at site i at depth t [t is the layer index). Then the transmission of weights can be represented via the stochastic 
equation, Eq. (|^) where it is assumed, for convenience, that the successive layers in the t direction are shifted by one 
lattice unit to the right. The injection term 1 represents the weight of a bead itself (assuming that all beads have the 
same weight unity) and represents the fraction of the weight that is transmitted from a given bead to its descendent 
in the next layer to the right. The only difference with the TM is that, in the q model, r^'s are independent and 
identically distributed random variables drawn from a uniform distribution over [0, 1]. Indeed one can study a general 
stochastic equation such as Eq. (|l|) where r^'s are independent and identically distributed random variables drawn 
from a general distribution /(r) over [0, 1] ||]. TM is a special case with /(r) = ^(J^.o + Similarly q-model with 

uniform distribution is another special case when /(r) = 1. 

For distributions of the form, /(r) = r"(l — r)"/B{n + l,n + 1) (where n is a positive integer and B(m,n) is 
the Beta function). Coppersmith et. al. showed |^ that the joint probability distribution of the normalized weight 
variables, Vi = rrii/t, factors in the limit t — s- oo, i.e., P(vi, V2, V3 . . .) = Yii Pi^^i) as t — > 00. The uniform distribution 
falls in this category as it corresponds to n = 0. Thus the correlations between normalized weights vanishes in the 
t — !■ 00 limit and the single point weight distribution was shown to be, P{v) — a(n)t;""^^ exp(— 2nu) for all v where 
a{n) is independent of f Experiments on bead packs measured the force distribution on the bottom layer of the 
pack (t —^ 00) and the results were found to be in agreement with the g-model results with n = 0, i.e., with uniform 
distribution, /(r) = 1 

While the spatial correlations between the normalized weights vanish in the t 00 limit, they are expected to 
be nonzero at finite depth t. We will show later that Claudin et. al. jl^ stated incorrectly that for the uniform 
distribution, the correlation is zero at any finite depth. Claudin et. al. also calculated |l^] the equal-depth correlation 
function in a continuum version of the g-model for generic distribution of the fractions and found a rather structure 
less correlation function. This is not surprising because they made the assumption that the beads are massless. Our 
exact calculation for the discrete g-model in this paper shows that for nonzero bead mass, the equal-depth correlation 
function has a very interesting scaling behavior characterized by a universal scaling function which is independent 
of initial conditions for short-ranged initial conditions. Besides, we also compute exactly the nontrivial temporal 
correlations between masses in the vertical direction. These temporal correlations have not been computed for the 
g-model before. 

C. New Results 

The new results that we obtain in this paper can be summarized as follows: 
(1) For the TM, we calculate exactly the equal-time mass-mass correlation function C{r,t) = {mo{t)mr{t)) — 
{mo{t)) {nirit)) between two spatial points separated by a distance r in all dimensions. We show that in the scaling 
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limit r oo, t ^ oo but keeping r/y/i fixed, C{r,t) = —r'G{r/\/t) where 7 = 2 for d < 2 and 7 = (3 — |) for 
d> 2 and the scahng function depends on d. In c? = 2, there are additional logarithmic corrections. In c? = 1, we also 
compute the scaling function G{y) explicitly. 

(2) Putting r = limit in the explicit expression for C(r, t) allows us to compute the on-site variance, {m?{t)) exactly. 
For large t, we find (w^) ~ t(4+'i)/2 in d < 2, (m^) ^ t^/\ogt in d = 2, and (m^) ^t^ md>i. For d = 1, this result 
was already derived by Takayasu et. al. but for d 7^ 1 it is a new result. This therefore proves rigorously, without 
any assumption of scaling, that the upper critical dimension (beyond which mean field exponents are correct) of the 
TM is dc = 2. 

(3) We also study exactly the normalized unequal time correlation function, Arit^r) = {Xo{t)Xr{t + 
T))/\/{^iit)){X?it + T-)) where Xr{t) = m,.(i) - {mr{t)), in all dimensions in the TM. The normalized autocor- 
relation function is obtained by putting r = r/2 in Ar(t,T). This is because each mass in the TM has a net drift 
velocity equal to 1/2 to the right. We show that the autocorrelation function, A^/2{t, r) ~ T~'^/'^h{T /t) in the scaling 
limit, r — > cx3, i — > 00 but keeping t /t finite. In d = 1, we derive the scaling function h[y) explicitly. 

(4) We also calculate exactly the correlations between forces at two different points (both equal depth and unequal 
depth) in the g-model |^ of force fluctuations in bead packs for arbitrary distributions of the fractions of weights. 
These correlations have so far not been measured experimentally. But if experiments can be performed in future, 
then our exact results will be useful for validation of the g-model. 



III. CORRELATIONS IN ONE DIMENSION 



In this section, we calculate exactly the spatio-temporal correlation function in the TM in one dimension. Even 
though the evolution equation for the single point probability distribution of mass P(m, t) involves the joint two point 
probability distribution function P(mi, m2, i), the evolution equation for the two point correlation involves only other 
two point correlation functions. This simplifying aspect makes the correlation function analytically tractable in the 
TM. 

The parallel dynamics of the TM in 1-d is represented by the stochastic equation (|^), namely, mi(i + 1) = mi(i)(l — 
Ti) + mi_i(t)ri_i -|- 1. If r.i = 1 at time t, the mass mi at site i jumps to its right neighbour while ri — Q indicates 
that it stays at site i. The hopping of the mass at site i to i -I- 1 is described by the first term while the second term 
accounts for the mass at i — 1 hopping onto site i. The last term 1 indicates the injection of unit mass from outside 
at every time t. Averaging Eq. (|l|) over all possible histories (starting from a zero mass initial configuration) we 
immediately get {m)(t) = t. 



A. Equal Time Correlations in One Dimension 



The evolution equation for the equal-time correlation function between two space points i and j can be written 
down by multiplying Eq. (0) by mj{t + 1) and then taking an average over all possible histories. Due to the 
translational invariance in an infinite lattice, this correlation function depends only on the difference |j — j|. Denoting 
X = i — j and using the translational invariance, we find the connected part of the correlation function at time t, 
Cx{t) = {mo{t)mx{t)) — t^, obeys the equation, 

C:c{t + 1) - \{C^+i + 2C, -f C,_i) -f ^(Co + t')(2(5,,o - - 4,-i). (3) 

In obtaining the above equation, we have used the fact that Vi and rrij at time t are independent of each other for all 
i and j. The Eq. (||) can be solved exactly for arbitrary initial condition by the generating function method. It turns 
out, however, that the solution at large time t becomes asymptotically independent of the initial condition as long 
as the initial condition is short ranged. Without any loss of generality, we therefore start from the simplest initial 
condition when the mass is at every site. Let F{q,u) — Yl^=iYl^o^x{t)<fu^- Multiplying Eq. (||) by q^u^ and 
summing over x and t, one can express F{q, u) in terms of Ci{u) = X^t^o C'i(t)u* and Gq{u) = X^t^o ^o(^)'^*- Ci{u) 
can further be expressed in terms of Co{u) from Eq. (|^) by putting x — Q, multiplying by u* and summing over t. 
Thus, finally wc get the following expression for F(q,u), 

q ^2(1 + u){l-q)- 2(1 - u)^Coiu)\ 
" (1 - m)3 [4g - «(1 g)2] ' 
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where C'o(u) is yet to be determined. We determine Co{u) by noting that F{q,u) has two poles qi^2 = {2 — u ± 
2\J\ — u)/u. For positive values of u, |gi| > 1 while \q2\ < 1. This would imply that for fixed time, Cx{t) will blow up 
exponentially as \q2\^ for large x. Since this can not happen, the numerator on the right hand side of Eq. (Q) must 
also vanish at q = (72 in order to cancel the pole. This immediately determines Co{u), 



A f ^ m(i + m)(i- yi-u) , . 

= — ir^^2 — ■ (5) 

Substituting Cq{u) in Eq. (^), the generating function F{q,u) is then fully determined, 

F{q,u) = — 7^ > , ,„ • 6 

Let us denote Cx{u) = X^t^o ^xit)u'^- The coefficient of q^ for a; > 1 can be easily pulled out from Eq. to yield 

~ -u{l + u) 

[\-uf + ^ ^ 

Note that it is evident from the above expression that Cxii) = for x >t. 

In order to derive the explicit expressions for Cx{t) for a; > 0, we need to invert the discrete Laplace transforms 
Cx{u). We first derive Co(i) explicitly by computing the coefficient of w* in the expression of 6*0 (u) in Eq. (^), 

2t{2t+l){At + l) {2t)\ 2 
4^5 

This therefore gives us an exact expression for the on-site mass variance Co(t) for all t. Taking the large t limit in 
Eq. (|) we get. 

We note that since Eq. (||) implies that (m^) ^ i(3-'^)/(2-'r)^ -^e get t = 4/3 in 1-dimcnsion. This therefore constitutes 
an alternate method to derive t = 4/3 in 1 dimension. 

In order to derive Cx{t) for x > 0, we need to calculate the coefficient of u* on the right hand side of Eq. (|^). For 
arbitrary x, this is somewhat hard. However it is easy to derive the asymptotic behavior of Cx{t) for large x and 
large t but keeping x/ y/t fixed. By taking u 1 limit in Eq. (|^) and after a few steps of algebra, we find that in this 
scaling limit 

Cx{t) = -t^G{x/Vt), (10) 

where the scaling function G{u) is universal, i.e., independent of the initial condition as long as the initial condition 
is short ranged and is given by 



fOC fOC fOC fOO 

G{y) = 32 dxi I dx2 / dx^ I dx4 erfc( 

J y J xi J X2 -^3:3 



0:4). (11) 

The complementary error function is defined as, erfc(y) = exp{—x^)dx. The above integrals can be done to 

derive an explicit expression for the scaling function. 



Giy) = l 



(3 + 12y2 + 4y4)erfc(y) - ^y(5 + 2y^)e-y" 



(12) 



We also performed a numerical simulation of the TM on a one dimensional lattice with periodic boundary condition. 
In Fig. 1, we show the scaling plot of the connected part of the correlation function. The data at different times, 
when scaled as in E q. (|lO| ) collapse onto a single scaling function which is in excellent agreement with the analytical 



result given by Eq. (|] 
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B. Temporal Correlations in One Dimension 



In this subsection we compute exactly the two time correlations in the TM in one dimension. Let us first define, 
Xx(t) = rrixit) — {mx(t)) — mx{t) — t. We then define the general two time correlation fimction as Dx{t,T) = 
{Xo{t)Xx{t + t)). It is also useful to define the normalized two time correlation function, 

Clearly then, {t,T) = (t,T)/ a/ Co(t)Co(i + r) where Co (t) is just the on-site mass variance whose exact expression 
is given by Eq. of the previous subsection and we have also used the translational invariance of the equal-time 
correlation function. Thus we just need to evaluate Dx{t,T) which can be done exactly as follows. 
From Eq. (|l|), it is easy to show that the function Dx{t,T) evolves as a function of r for fixed t as, 

D,{t, r + {Dx{t, t) + r)) , (14) 

starting from the initial condition D,j.{t,Q) ~ Cx{t), where Cx{t) is the equal-time correlation function computed 
already in the previous subsection. Let H{k,t,T) — J^'^^-oo -^x{t,T)e^'^^ . Note that here we used the x summation 
from — oo to oo as opposed to to oo. This is because Dx{t,T) is not equal to D^a:{t,T) for r > 0. They become 
equal only for r = due to translational invariance. From Eq. ( p^ ) we get 

H{k,t,T) ^ H{k,t,0) (-^j . (15) 

Note that H{k,t,0) = E^-oo Cx{t)e''''' and C^.x{t) = C^(t) as translational invariance holds for equal-time correla- 
tion function. 

By inverting the Fourier transform in Eq. (|l5|), we get a simple expression for Dx{t,T) in terms of the equal-time 
correlation fimctions, 

DAt,T) = ^ j^J dkH{k,t,T)e-'''- 

= (16) 

In order to calculate the auto-correlation function, we note that in the TM, the masses have a net drift velocity, 
V = 1/2 towards the right. This is because of the definition of the model: in one time step, a mass either stays at 
its own site with probability 1/2, or hops to the neighbour on the right with probability 1/2. Thus, to calculate 
the proper auto-correlation function, one has to compute it in the moving frame which is shifting towards right with 
uniform velocity 1/2. Hence the correct auto-correlation function would be given by 0^/2(1, t). Putting x = r/2 in 
Eq. ( p^ ) and taking the transform, Dt-/2{u,t) = X^t^o ^t"/2(^i ''")"*' S^^ 



Dt/2{u,t) = ^ 



^ ' ' m=T/2+l ^ ' 



(17) 



where we have used the symmetry Cx(u) = C-x{u). Using the exact expressions for Cx{u) from Eqs. (|^) and (^ of 
the previous subsection in Eq. (p^ and after some steps of algebra, we get 



r/2-1 



n r., - 1 ^ ^ "(1 + ") u{l + u) u{l + u) ^ (2i\^uy 



1=0 

where we have used the combinatorial identity |E(i, 



,,_ (l + fc)" fc-Zy n \ (l-fc)(l-fc)"-i "g;V2»^ ( k ^ 
3) " 2 2 \nl2) 2 ^ \^)\{l^kY) ' ^^'^ 



j=n/2+l ^ ' \ ' / 
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In order to analyze Eq. (18), we first put u = 1 — s and note that the equation allows a scaling limit when s 
and T oo but the product sr remains fixed. In terms of time variables, this scaling limit corresponds to r ^ oo, 
i — > oo but keeping the ratio r/t fixed. In this limit, we find 



-Dr/2(S,' 



1 



where the scaling function is given by 

9i{y) 



(20) 



(21) 



In terms of time variables r and t in the scaling limit, t — )■ cxd and t ^ oo but keeping t /t fixed, we get the following 
expression by inverting the Laplace transform in Eq. (po|). 



16V2 ft T 



where 



/.r(.)^^^ 6sin- 



1 \ Vyjw + '-^y) 

^/TT^ x/2(l + y/2)2 



(22) 



(23) 



Using the above result and the exact large t behavior of Co (t) from Eq. (^ in Eq. (jl^) , we finally obtain the scaling 
behavior of the normalized autocorrelation function in the scaling limit mentioned above. 



1 



where the scaling function is given by. 



Hy) 



TT (l + y)5/4 



l-^(l + y/2)^ 6sin-i 



1 \ yyiw + Sy) 

yTTW2/ V2(l + y/2)2 



(24) 



(25) 



The function h{y) ^ y as y — > and h(y) « y if^J/ as y — )■ oo. Thus, for large y, the scaling function decays 
as a power law with an exponent 9/4. 

We remark that the above scaling behavior holds only in the limit when r — > cx), i ^ oo but the ratio r/t is 
kept fixed. In other limits, it is also possible to investigate the detailed behavior of the autocorrelation function by 
analyzing the exact equation (18). 



IV. CORRELATIONS IN ARBITRARY DIMENSIONS 



In this section we study the two-point spatio-temporal correlations in the TM in an arbitrary spatial dimension d. 
As mentioned in section II, it turns out that for general d, equations for the correlations simplify considerably for the 
continuous time version of the TM. In this version, every mass hops with rate p to each of its d nearest neighbors in 
the positive direction, and aggregates with the mass present at the hopped site. In addition, injection of unit mass 
occurs at every lattice site with rate 1. 

The evolution of the mass m{xi,X2 ■ ■ ■ Xd, t) in a small time interval At can be represented by the equation, 

d d 

m{{xi},t + At) = ^rj"m(xi, . . . ,Xj - 1, . . . ,Xd,t) + {1 - ^ rj + )TO({a;J, i) + /({a;J,i), (26) 

where r^^'s are independent and identically distributed variables, with distribution /(r) — pAtSr^i + (1 ~ pAt)Sr,o 
and indicate the hopping events of the particles. The random variable I{{xi},t) denotes the event of injection and is 
drawn from the distribution, P{I) — AtSj^i + (1 — At)6i^o independent of the r^^'^'s. 



7 



A. Equal-Time Correlations 



From the equation ( |2q ) of evolution of the masses, one can easily write down the evolution equation for the two-point 
correlation function in continuous time. Multiplying Eq. (|2^) at two different space points and neglecting terms of 
order O(At^) and higher, we find that the two-point correlation function, C{{x}, t) = (m(0, . . . , 0, t)m{xi, . . . , Xd, t)) — 
t'^, evolves as 

d 

—C{{x})^-2pdC{{x})+p^ X! C{xi,. . . ,Xj +m,. . . ,Xd) + S.j,^fi . . .S.j,^fi 

j—l 7n—±l 

d 

+ P (C({0}) + 1^) {2dS^,^o ■ ■ ■ 5^afi - 

J = l 



(27) 



where we have suppressed the t dependence of C({x}, t) for notational convenience and also used translational invari- 

C({x}) exp(i ^^^j^ kjXj) then evolves as. 



ance of C{{x},t). The Fourier transform G({fc},f) — X]{^}=-oo 
d 



dt 



j=i- 



G{{k}) = 2p (G({fc}) - cm) -t')\-d + J2 cos(fc, 
Taking Laplace transform with respect to t, we get 

Fi{k},s)^ 



(28) 



s2 + 2p(s3go(s) + 2) 


^c^-Ej=iCos(fcj-)) 


s3 


s + 2j3(^d-X;^^iC0s(fcj)) 





(29) 



where F{{k}, s) = G{{k}, t)e ^*dt and go{s) = C(0, . . . , 0, t)e ^*dt which is yet to be determined. We deter- 
mine go{s) by noting that go{s) — lo^ ■ ■ ■ lo^ -^({^Ij s)dki . . . dkd- Integrating Eq. (|9|) with respect to the fc^'s 
we get the following expression for ga{s), 



5o(s) = ^(s'-2s+-^) 



where 



/(.) = 



(2^)^ Jo 



2-K 



dk\ 



dkd 



(30) 



(31) 



s + 2p\d-Y.j=\^^'^{^i) , 

The small s behavior of /(s) can be easily evaluated by analyzing the integral in Eq. (pT|). We find that as s — > 0, 

/(s) - s~(i~5) for d < 2 
log(s) for d = 2 

^ constant for d > 2. (32) 

Substituting /(s) in the expression for go{s) in Eq. ( ^0| ) and inverting the Laplace transform, we find that the on-site 
variance for large t behaves exactly as, 

C(0, ...,0,i) for d<2 

t3 



for d = 2 



log(t) 

for d>2. 



(33) 



Note that the results in Eq. (^) are exact results for large t and does not assume any scaling behavior. This result 
clearly proves rigorously that the upper critical dimension of TM is dc = 2. 

In addition, if we assume that the on-site mass distribution scales as, P{m,t) ^ m^'^ f{mt~^) for large m and large 
i, we get the following results for r and 5. The first moment, (m) ~ t gives (5 = 1/(2 — r) [|^. The second moment 
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scales as, (m^) t^^~'^'^^. Using the exact results for variance from Eq. (|33|), we get, r = 2(d + l)/{d + 2) for d < 2, 
T = 3/2 ior d > 2 and T = 3/2ind = 2 with additional logarithmic corrections. These results for r are in agreement 
with the results obtained by Swift et. al. |Q by using a more indirect mapping to the age distribution of particles in 
a related reaction-diffusion process and also assuming scaling behavior. 

With go{s) completely determined, we therefore have an exact expression in Eq. ( p9| ) for F{{k},s), the joint 
Laplace-Fourier transform of the full correlation function C{{x},t) in arbitrary dimensions. For arbitrary d, it is 
complicated to invert this transform to obtain an exact expression for C'{{x},t). However, by analyzing the small s 
and small k behavior of F{{k},s), it is easy to see that for large x and large t, but keeping fixed, C{{x},t) 

satisfies a scaling behavior, C{{x}, t) ^ f^Gil^}) with 7 = 2 for d = 1 and 7 = (3 — |) for c? > 2 and the scaHng 
function G{y) depends explicitly on d. 

For d = 2, there is additional logarithmic correction and the scaling breaks down. In this case, the exact expression 
for /(s) is given by, 

K{4p/{s + 'lp), (34) 



7r(s + 4p) 

where K is a. complete elliptic integral pG]. This gives an explicit expression for go{s), 



7r(s -t- 4p) 12 



Substituting the expression for 50 (s) in Eq. (|29|) we get, 

2 7r(s-t-4p) 
^ ^'^'^7^' s^K{Ap/{s + Ap)) [s + 2p{2 - cos(fci) - cos(fc2))] ' 

After some straightforward algebra, it turns out that the large distance behavior of C{x,y,s) = C{x,y,t)e^''*dt 
is given by, 

r<f \ -(s + 4p) f ry/s \ 



where r = \/x'^+y'^ and Kq is the modified Bessel function pO[ . In order to get an explicit expression of C{x,y,t) 
for large r and t, one needs to invert the Laplace transform given by Eq. (|37|). But it is obvious from this expression 
that C{x,y,t) will no longer have a nice scaling form in the large distance and long time limit as in one dimension 



(see Eq. (10)) due to the appearance of logarithms in the asymptotic behavior of the functions K{x) and Ko{x). This 



violation of scaling due to logarithms is again expected since 2 is the upper critical dimension of the TM. 



B. Temporal Correlations 

We can write down the equations for the time evolution of the temporal correlation function in a manner very similar 
to that in d = 1. We define the connected correlation function as I?a;j^...^a;^(t, r) = {mQ^,,,fi(t)mx-^^,,,^,j;^(t + T)) — t{t + T). 
From the evolution equation of the masses, it is easy to show that the function Dx{t, r) evolves as a function of t for 
fixed t as, 

^i?{,}(i,T) =p |^^(7?,,.....,^._i,....,,(i,r)-d7?{,}(i,r)j , (38) 
starting from the initial condition D^^^{t^()) = C^x^{t) where C^x}{i) is the equal-time correlation function. Let 

.{x}-- " 



H{{k}, t, r) = J2t^=~oa D{x}{t, t) exp(V ^0 ^kjXj). From Eq. (||) we immediately get 



t, t) = H{{k}, t, 0) exp I pT - 1) I , (39) 



where H({k}, i, 0) is the Fourier transform of the equal-time correlation function. 
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We invert Eq. ( |39D to get the temporal correlations in terms of the equal-time correlation functions, 

d 



„27r d oo d d 

/ d''/cexp(-i^%Xj) ^ C({a;'},i)exp(i^fcjx^)exp(pr^e*'=0. (40) 

"'0 r..,i .-I 



(2^)^ 7o 



j = l {x'} = -oo j=l j=l 



As in d = 1, there is a net drift velocity, pr, in each forward direction. Hence the correct autocorrelation function 
is given by £'{pr}(t, t). Putting xj — pr in Eq. (|40|), and simplifying, we get 



D{pr}{t,T)^ V C(pT-mi,...,pT-m,)e-'^P^^^^— . (41) 



It can then be shown that in the scaling limit, t — > oo, t oo but keeping r/t fixed, the normalized autocorrelation 
function allows for a scaling solution as in d = 1 , 

^{Pr}it^^)-^H^)- (42) 

We do not present the explicit form of the scaling function here. It can be shown that h{y) — > const, as y — s- 0. Thus 
in arbitrary dimensions, the asymptotic decay of the normalized auto-correlation function is given by, t^'^/'^ for large 
T (after taking the large t limit). 



V. CORRELATIONS IN THE q-MODEL OF FORCE FLUCTUATIONS 

The g-model of force fluctuations has been defined in section (II. B). In this model, the variables evolve in time 
via the stochastic equation, mi{t -f 1) = (1 — ri)mi(t) -\- ri_imi-i{t) + 1, where the random variables r^'s are drawn 
independently from a arbitrary distribution f(r) over the support [0, 1]. Experimental results for the force distribution 
in real bead packs were found to be described accurately by the q- model with uniform distribution, f{r)~l [p^ . 

In this section, we calculate exactly the two point correlations between to^'s for the generalized g- model, i.e., for 
any arbitrary distribution /(r). It turns out that the two point correlations are characterized by two parameters, 
/ii = rf{r)dr and ^2 = Jq r^f{r)dr with > ^2- In the /ii > ^2 plane, there are two types of asymptotic 
behaviors depending on whether /ii = 112 or fj,i > ^2- At every point on the line fJ.i — fJ-2, the correlation function 
has the same universal asymptotic behaviour independent of initial conditions as long as they are short ranged. The 
special case of the TM with fii = fj.2 = 1/2 falls within this class. On the other hand, for all points in the plane 
fJ'i > ^2, the correlation function has once again the same universal asymptotic behavior regardless of the actual 
values of /xi and fi2 but this behavior is different from the behavior on the line fii — fj.2- The q- model with uniform 
distribution corresponds to the point /ii = 1/2, /i2 = 1/3 and therefore falls in the second category. 

Since the steps of the calculation follow closely that of the TM, we will skip most of the details and present only 
the final results. 



A. Equal-Time Correlations in One Dimension 

Starting from the stochastic evolution equation (Eq.(|l|)) of masses in d = 1, it is straightforward to write down 
the evolution equation for the equal time correlation function, Cxit) = {mo{t)'inx{t)) — t^. We find that for general 
distribution /(r), the equal-time correlations evolve as 

a(t -f 1) = (mi - Mi')(a.+i + c,-i) + (1 - 2^1 + 2m?)c, 

+ {^l2 - Mi')(C^o + i')(24,o - 4,1 - 4,-1), (43) 

where /ii = rf{r)dr and ^2 — r^f{r)dr. Note that for the TM, /ii = 1/2 and /i2 = 1/2 and then Eq.(^) reduces 
exactly to Eq. (||) studied in section (III. A). For the imiform distribution, /(r) = 1, one gets, /ii = 1/2 and fi2 — 1/3. 
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We can solve the Eq.(^3|) for arbitrary initial condition and arbitrary parameters /zi and fjL2- It turns out that the 
dependence on the initial condition drops out for asymptotically large t. Following the same steps as in the TM in 
section-Ill. A, we find that the Laplace transforms, Cx{u) = X^t^o C^x{t)u*, are given exactly by 



^ , . u{l + u)[^l + (4a - 1)^ - VT^] 

L^0(.?ij = J-. 7^ ]—{ , (,44 j 

(1 - uYgi{u) 

\ 2x 

^ _ n., _ . /I _ \ 

(45) 



(1 - u)2gi{u) 




where gi{u) ~ [^-j-^-y/l + (4a — l)u + y^l — u] and a = /^i — /ii^ and 6 = /i2 ^ Mi^- It is clear from the above expressions 
that the asymptotic behaviors for large t (corresponding to u — s- 1) depend on whether a = b, i.e., fj,i = 112 or a > 6, 
i.e., /zi > ^2- 

For /ii = /i,2, we find, after inverting the Laplace transforms, that for large t, 

C,(t)»^t^/\ (46) 
Cx{t) « -fG (^=] , x,t» 1, (47) 



/4at 

where the universal scaling function G{y) is given by Eq. ( p^ as calculated for the TM. 
For /ii > /i2, on the other hand, we find for large 

Co(0 ~ -^t\ t » 1, (48) 
a — b 

Cx{t) « -^j^—tic, (^) , x,t» 1, (49) 



y/a{a~b) \^/4M 
where the universal scaling function Gi{y) is given by 

Gi(y) = \ (^-^—^^^ _ y[i + 2y')cMy)^ . (50) 

Note that the uniform distribution corresponds to /ii — 1/2 and /i2 — 1/3, i.e., a — 1/4 and b — 1/12. In Fig. 2, we 
compare the numerically obtained scaling plots for the uniform distribution with the exact scaling function given by 
Eq. (H). 

If we define the scaled weight v ~ m/t as in Ref. ||], then the connected part of the two-point correlations 
{v{Q^t)v{x,t)) — 1 ~ — i^^/^Gi(a;/\/4at) for large x and t. Clearly as t ^ 00, the scaled weights w's get completely 
uncorrclatcd but for finite depth (or time) t, there is a nonzero anti-correlation specified by the scaling function Gi{u) 
which might be possible to measure experimentally. We also point out that the statement of Claudin et. al. ||l9| that 
for the uniform case, the correlation is zero at any altitude t is clearly incorrect. 

We note that from Eq. (^3|), one can easily derive the evolution equation for the correlation function in the 
continuum space and time. For a proper continuum limit in time, we need to assume that both a — > aAt and b — > bAt 
are of order At. Defining the Fourier transform, G{k,t) — j'^^G{x,t)e'^^^ dx, one finds from Eq. ( ^3|) that for small 
k, the correlation function evolves via the equation. 



{dt + ak^)G{k, t) ^bk^(^J ^G(k', t) 



+ t^j. (51) 

Note that this equation above is identical to the one derived by Claudin ct. al. |l^ except for the additional k^t^ 
term on the right hand side of Eq. (|5l|). The origin of this additional term can be traced back to the fact that in our 
case, the bead mass is nonzero (equal to 1) as opposed to the zero mass case considered in reference [|l^. As seen 
from our analysis that for nonzero mass, the correlation function has a much more nontrivial and universal structure 
as opposed to the rather structure less correlations found in the zero mass case in [ p^ . 

Rescaling time by the factor a, one sees that Eq.(|5l|) is parameterized by the single variable X — b/a. As in the 
discrete case, there are two possible asymptotic behaviors depending on the value of A. For A = 1, one finds the 
Takayasu type behavior and a completely different asymptotic behavior emerges for all A < 1. 
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B. Temporal Correlations in One Dimension 



The two time correlations can similarly be computed for the g- model for any arbitrary distribution /(r). We present 
here only the explicit results for the uniform distribution in d = 1. 

For two time correlations we use the same notation as in the TM (see section (III.B)). For the uniform distribution, 
/(r) = 1, we find that in the scaling limit, t —^ oo and i — > oo but keeping r/t fixed, the normalized auto correlation 
function, as defined in Eq. (|l3|), has the scaling behavior, 



(52) 



where the scaling function is given by 

h2iy) - 



The scaling function h2{y) 



as y 



7r3(l + y) 

as in the TM. For large y also, h2{y) decays as a power law, h2{y) 



(53) 



2 as in the TM but with a different exponent 2 than for the TM exponent 9/4. 



n{k},s)^ -^Z."'T' ^ (54) 



C. Results for Arbitrary Dimensions 

Following the similar line of arguments as in the TM, the evolution equation for the two-point correlation function 
can be derived for the generalized g-model in arbitrary spatial dimension with discrete space and time. The discrete 
equations are rather complicated but the asymptotic properties can be easily derived by taking continuum space 
and time limit. In the continuum limit, it turns out that as in d = 1, the Fourier transform in any dimension, 
G({fci}, t) = C{{xi}, t) exp(i ^ kiXi)d'^x evolves by the simple equation ( ^l|) parameterized by the ratio \ — b/a 
once time is rescaled by a. The parameter A < 1 as in the discrete case. As in the TM, taking the Laplace transform 
of Eq. (Isi) ) with respect to t, we find 

Afc2[ 2 + ,9o(g)g 3i 
5(fc2 + s) 

where F({fc}, s) — G{{k}, t)e~'^*dt and go{s) — C(0, . . . , 0, t)e^'^^dt which is yet to be determined. As in the 
TM, go{s) is determined by integrating Eq. ( p4| ) with respect to k. Note that the upper cut-off of each ki integration 
is now set by 27r/A where A is the lattice constant. We find 

2\[l~ sl{s)] 
s-'[l — A -f As7(s)J 

where I{s) = J d'^k/ (fc^ -|- s). The small s behaviour of I{s) is same as given by Eq. (|3^). 

From Eq. ( |55| ) it is clear that there are two different small s behaviors of go{s) depending on whether A = 1 or A < 1. 
The TM corresponds to A = 1 and its asymptotic behaviors have already been discussed in detail in Section- (IV). 
Here we focus on A < 1. In that case, using the small s behavior of /(s) in Eq. (|5^), it is clear that 50 (s) ^ 2X/ in 
any dimension, indicating that C(0, 0, . . . , 0, f) ^ for large t. Thus for A < 1, in contrary to the A = 1 case (TM), 
there is no critical dimension separating different asymptotic growth of C(0, 0, , . . . , 0, i). 

Substituting this expression for gQ{s) in Eq. (p4), we find that for small s 



0\ 7,2 

^(«'^)-(r3WT^- ^''^ 

It is then not difficult to derive the asymptotic properties of the correlation function in real space and time. We find, 

C{{Q},t) r-^e, t»l, (57) 
C({x}, t) « -t(2-|)g^ ^ ^ ^^^ (53^ 

where Gi (y) is the dimension dependent scaling function. 

The most important result of this subsection is that while for A = 1, there is an upper critical dimension dc — 2 
separating different asymptotic growth of the on site variance, there is no such critical dimension for A < 1 (which 
includes the uniform distribution of fractions for the g-model) . 
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VI. SUMMARY AND CONCLUSION 



In this paper we have computed exactly both the equal-time as well as the unequal-time two point correlations in 
the Takayasu model of mass aggregation and injection in all dimensions. We have identified different scaling limits 
and obtained the scaling functions explicitly in d — 1. Our exact results for the on-site mass variance prove rigorously, 
without any assumption of scaling, that the upper critical dimension of the Takayasu model is dc = 2. 

We have also extended our technique to compute exactly the correlations in a larger class of aggregation models 
with injection. This generalized model includes, as special cases, the Takayasu model and also the g-model of force 
fluctuations in granular materials. We have shown that the correlation functions in this generalized model is parame- 
terized by two variables (/ii, IJ.2) which are respectively the first and the second moment of the distribution f{r) of the 
fractions r^'s. We have shown that in the two-dimensional parameter space (fii, with fj,i > /i2, there are two types 
of asymptotic behaviors of the correlation function depending on whether fii = ^2 or fix > ^2- For generic points in 
the region fii > fi2 which includes the uniform distribution represented by (1/2, 1/3), the correlations have similar 
asymptotic behaviors which is different from that on the line /ii = 112 which includes the Takayasu model. Besides, 
for fii > ^2, there is no upper critical dimension in contrast to the case ni = ^2 where the upper critical dimension is 
dc~2. We have presented explicit forms of scaling functions for both the Takayasu line as well as the experimentally 
relevant uniform distribution case. These exact results will be useful for comparison with possible future experimental 
results on correlations in bead packs. 

In this paper we have calculated exactly various time-dependent correlations between forces in bead packs in the 
context of the simple scalar g-model. There have been various generalizations of this scalar model to include the 
tensorial nature of the forces [|l9|,^ and also to non-cohesive granular materials It would be interesting to see 
if our method can be extended to calculate the correlations in these generalized models. 

We thank D. Dhar and M. Barma for useful discussions. 
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0.5 1 1.5 2 2.5 

FIG. 1. The figure shows the scaling plots of the equal-time correlation function, Cx{t) — {m{0,t)m{x,t)) —t^ obtained from 
numerical simulation of the TM on a one dimensional lattice of 1000 sites. The data for five different times collapse onto a 
single scaling curve which matches very well with the analytical scahng function given by Eq. (|l2). 
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xj \/t 

FIG. 2. The figure shows the scaling plots of the equal-time correlation function, Cx{t) — {m{0,t)m{x,t)) — t^ obtained from 
numerical simulation of the g-model on a one dimensional lattice of 1000 sites. The data for five different times collapse onto 
a single scaling curve which matches very well with the analytical scaling function given by Eq. (H). 
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